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From direct A-body simulations we find that the dynamical evolution of star clusters is strongly influenced by the 
Roche volume filling factor. We present a parameter study of the dissolution of open star clusters with different Roche 
volume filling factors and different particle numbers. We study both Roche volume underfilling and overfilling models and 
compare with the Roche volume filling case. We find that in the Roche volume overfilling limit of our simulations two- 
body relaxation is no longer the dominant dissolution mechanism but the changing cluster potential. We call this mechnism 
“mass-loss driven dissolution” in contrast to “two-body relaxation driven dissolution” which occurs in the Roche volume 
underfilling regime. We have measured scaling exponents of the dissolution time with the two-body relaxation time. In 
this experimental study we find a decreasing scaling exponent with increasing Roche volume filling factor. The evolution 
of the escaper number in the Roche volume overfilling limit can be described by a log-logistic differential equation. We 
report the finding of a resonance condition which may play a role for the evolution of star clusters and may be calibrated 
by the main periodic orbit in the large island of retrograde quasiperiodic orbits in the Poincare surfaces of section. We also 
report on the existence of a stability curve which may be of relevance with respect to the structure of star clusters. 
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1 Introduction 

While open clusters (OCs) dissolve, they are populating the 
disc with field stars by continuous mass loss in the tidal field 
of the Milky Way. The fraction of field stars which belonged 
once to a star cluster is estimated to be 10% or less (Miller 
& Scalo 1978; Wielen 1971), up to 40% (Roser et al. 2010) 
or 100% (Maschberger & Kroupa 2007). 

The feeding of the field population by OCs is deter¬ 
mined by the initial cluster mass function (ICMF), the clus¬ 
ter formation rate (CFR) and the mass loss rates of clus¬ 
ters corrected for stellar evolution. The current ICMF can 
be directly determined by the young clusters of an unbiased 
cluster sample in the solar neighbourhood (Kharchenko 
et al. 2013) or in nearby galaxies (Boutloukos & Tamers 
2003). But for disentangeling the CFR and the cluster 
mass loss rates in the cluster distribution function N{M, t ) 
with mass M and age r additional information is needed. 
Tamers, Gieles & Portegies Zwart (2005a) used the dissolu¬ 
tion timescales of Baumgardt & Makino (2003) to analyse 
N{M,t) in the TMC, M33, M51 and in the solar neigh¬ 
bourhood. In Tamers et al. (2005b) a more detailed analysis 
of the solar neighbourhood has shown that the present day 
CFR provides 30% of the star formation rate as determined 
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by Just & JahreiB (2010). It was also shown in Tamers et 
al. (2005b) that the dissolution timescales of OCs in the so¬ 
lar neighbourhood are significantly smaller than those pre¬ 
dicted by Baumgardt & Makino (2003). 

A major issue in the determination of cluster mass loss 
and lifetimes is the wide range of possible initial conditions. 
OCs form in giant molecular clouds and after a few Myr 
the volume is cleared from the interstellar medium. In this 
phase a large fraction of embedded clusters may become un¬ 
bound and dissolve quickly (sometimes called “infant mor¬ 
tality”). Since the dynamical state of the stellar component 
is strongly affected by the loss of the cloud potential and 
because the Galactic tidal field is not dominating the clus¬ 
ter formation process, the initial conditions of the isolated 
young clusters are not well constrained. The structure of the 
cluster may be characterized by three aspects, namely the 
dynamical state, the stellar content, and the density profile. 
After the gas removal, the OC may be supervirial and ex¬ 
pand in a violent relaxation phase (e.g. Parmentier & Baum¬ 
gardt 2012). This phase lasts for a few crossing times and 
a compact core may survive, which then builds the start¬ 
ing point of the longterm evolution of the cluster. Therefore 
most investigations of the dynamical evolution of star clus¬ 
ters start in dynamical equilibrium and then add the tidal 
field. An initial mass segregation may alter the dynamical 
evolution of the cluster. Mass segregation is observed in 
some young, massive and compact clusters (e.g. Habibi et 
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al. 2013; Pang et al. 2013). Additionally, mass loss by stellar 
evolution depends strongly on the shape of the adopted ini¬ 
tial mass function (IMF). A pioneering “survey” of the dis¬ 
solution of stellar clusters is the work of Fukushige & Heg- 
gie (1995). They investigated a parameter space of different 
concentrations and slopes of the IMF and found that less 
concentrated clusters and/or those containing more massive 
stars are disrupted sooner. 

The density distribution of a cluster may be character¬ 
ized by the core radius Tc, the half-mass radius r/j, and 
the cutoff radius r*, where the density drops to zero. The 
general shape of the density profiles can be measured by 
the Lagrange radii ri% containing i% of the cluster mass 
(rh = r’ 50 % and rt = rioo%)- In the widely used King mod¬ 
els (lowered isothermal spheres, King 1966) the concentra¬ 
tion c = log(rt/rc) (or equivalently the depth of the poten¬ 
tial well Wo) is a free parameter, which can be set to any 
positive number (see Binney & Tremaine 2008, for more 
details). In contrast, the ratio of cutoff to half-mass radius 
varies only by a factor of ~ 3.3 with the minimum at low 
concentration (r^/r/j = 3 for Wq = 1) and a maximum 
at Wo ~ 8 (rt/vh = 9.1), which turns out to be a serious 
restriction for setting initial conditions of extended clusters. 

The strength of the Galactic tidal field can be quantified 
by the Jacobi radius rj, which is the distance of the La¬ 
grange points Li and L 2 , the saddle points of the effective 
potential, to the cluster centre. The Jacobi radius for circular 
orbits is given by (Just et al. 2009; Kiipper et al. 2008) 


GM,i{rj) 


where G, Mci{rj), ftc and j3c = kc/^C the grav¬ 
itational constant, the cluster mass inside rj, the circu¬ 
lar frequency and the ratio of epicyclic over circular fre¬ 
quency, respectively. For clusters on a circular orbit (in an 
axi-symmetric Galactic potential) the Jacobi energy Ej is 
a constant of motion and all stars inside the Roche volume 
given by rj with Ej < i?j,crit, where i?j,crit is the effec¬ 
tive potential at the Lagrange points Li and L 2 , cannot es¬ 
cape. But there is no general criterion for a bound system. 
Fellhauer & Heggie (2005) have shown that unbound, low 
density systems can survive for more than a Gyr in the tidal 
field. In general potential escapers with Ej > £^j,crit may 
stay a long time in the vicinity of the cluster, before they es¬ 
cape to the tidal tails or return to the cluster. Ross, Mennim 
& Heggie (1997) derived a criterion for escape dependent 
on the offset of the guiding radius and the epicyclic radius 
of the star orbit. Applying this criterion to a flat rotation 
curve, the closest point of ‘safe’ escape is at a distance of 
~ 2.6r,/. They have also shown that stars with arbitrarily 
large epicylic motion may return to the cluster. Therefore it 
is appropriate to count all stars inside 3rj to be bound as a 
simple criterion. 

For this study, we will use the 100% Roche volume fill¬ 
ing factor A = rioo%/r’j to set the initial size of the clus¬ 
ter with respect to the tidal field. For practical reasons (see 


Section 2) we will use A' = rgg%/rj to scale the initial 
size of the cluster and determine A analytically. Since the 
outer shells of the cluster contain only a small fraction of 
the cluster mass, the half-mass Roche volume filling factor 
A = Th/rj is a more robust measure of the impact of the 
tidal field on the cluster evolution. Our goal is the analysis 
of the evolution of Roche volume overfilling star clusters. 
There is a huge number of publications on numerical sim¬ 
ulations concerning the dissolution of Roche volume filling 
or underfilling star clusters in tidal fields. We can refer only 
to a selected subset representing the main results relevant 
for our purpose. 

Engle (1999) examined for the first time the evolution of 
Roche volume underfilling star cluster models and found a 
relaxation driven expansion phase until the previously un¬ 
derfilling clusters filled the Roche volume. Engle (1999) 
also presented for the first time a plot of lifetime vs. recip¬ 
rocal Roche volume filling factor from direct A^-body sim¬ 
ulations and found increasing lifetime for decreasing filling 
factor. 

Eukushige & Heggie (2000) calculated the Jacobi en¬ 
ergy dependence of the time scale of escape from a star 
cluster in a tidal field using a theoretical result of MacKay 
(1990). 

Baumgardt (2001) used the calculation by Eukushige 
& Heggie (2000) and obtained the scaling of the disso¬ 
lution time with the two-body relaxation time for poten¬ 
tial well filling clusters in a tidal field. Eor the latter case 
he postulated a steady state equilibrium between escape 
and backscattering of potential escapers into the poten¬ 
tial well and found that the half-mass time scales with 
[A^/ 10 ( 7 ^)]^^"^ in that equilibrium, where 7 is the factor 
in the Coulomb logarithm occuring in the two-body relax¬ 
ation time given by Eqn. (3) below. The potential escapers 
are stars which have been scattered above the critical Jacobi 
energy but which have not yet left the star cluster region. 

Baumgardt & Makino (2003) presented the results of 
a large parameter study of the evolution of multi-mass star 
clusters in external tidal fields (i.e. a logarithmic halo). They 
used different particle numbers, orbital eccentricities and 
density profiles for star clusters, tending more in the direc¬ 
tion of globular clusters rather than towards the regime of 
OCs. 

Tanikawa & Eukushige (2005) simulated a comprehen¬ 
sive set of equal-mass star cluster models with different 
Roche volume filling factors including for the first time 
Roche volume overfilling clusters. They cover a wide range 
of particle numbers reaching the globular cluster regime and 
quantify how the dependence of the mass loss time scale on 
the two-body relaxation time scale depends sensitively on 
the strength of the tidal field as imposed by the Roche vol¬ 
ume filling factor. 

If two-body relaxation drives the evolution, the term 
proportional to [A^/ 10 ( 7 ^)]^^"* in the half-mass time ac¬ 
cording to the theory in Baumgardt (2001) can be well 
approximated by BN'^ with p « 0.6 (Earners, Gieles 
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& Portegies Zwart 2005a). Lamers, Gieles & Portegies 
Zwart (2005a), Lamers et al. (2005b), Gieles & Baum- 
gardt (2008) and Lamers, Baumgardt & Gieles (2010) find 
that the dissolution time (e.g. half-mass time) scales with 
for the Roche volume filling case, where 
Mci is the initial cluster mass and the exponent depends on 
the parameter ILq of the King (1966) initial model. 

This particularly means that the equilibrium postulated 
by Baumgardt (2001) may not be realized. Initially this can 
be in fact true, as Baumgardt himself notes. The reason is 
that the potential escaper regime may be initially overpopu¬ 
lated. The dependence on Wq noted above is a hint that this 
scenario plays a role. Baumgardt (2001) did not state the 
linear stability analysis of the equilibrium which he postu¬ 
lated in his 2001 paper. The proof that the postulated steady 
state equilibrium is attracting typical non-equilibrium ini¬ 
tial states and the quantification of the attraction strength 
and time scale are still open issues. 

In the present study, we will present numerical evidence 
for the fact that open star clusters in the Roche volume over¬ 
filling regime dissolve mainly due to the changing cluster 
potential and the shear forces of the differentially rotating 
galactic disk. We call this mechanism “mass-loss driven dis¬ 
solution” in contrast to the “two-body relaxation driven dis¬ 
solution” which occurs from the Roche volume underfilling 
regime up to the Roche volume filling case (see also White- 
head et al. 2013, based on simpler models). 

We concentrate in the present study on the properties 
and evolution of classical OCs which already left their par¬ 
ent molecular cloud. Therefore the formation process and 
the early phase of gas expulsion of OCs is beyond the scope 
of the project. On the other hand, mass loss of the OCs 
due to stellar evolution (supernovae, stellar winds, plane¬ 
tary nebulae) will be taken into account. We systematically 
study the influence of the Roche volume filling factor which 
we have chosen as the main free parameter. In this sense, 
the present study aims to extend the study by Baumgardt 
& Makino (2003); Engle (1999); Tanikawa & Fukushige 
(2005) into the overfilling regime. In particular, we aim at 
obtainig scaling exponents of the dissolution time with the 
two-body relaxation time scale and at deriving a fitting for¬ 
mula for the dissolution time (half-number time). 

This paper is organized as follows; In Section 2 we 
shortly explain the method of direct A^-body simulations 
in an analytic background potential and the programs 
NBODY6tid and <y9-GRAPEH-GPU. In Section 3 we discuss 
the parameter space. Section 4 contains the results and Sec¬ 
tion 5 the conclusions. 

2 Method 

The dynamical evolution of OCs is calculated as an N- 
body problem in an analytic background potential of the 
Milky Way. For the background Milky Way potential, we 
use the same model as in Ernst et al. (2010, 2011); Just et 
al. (2009); Kharchenko et al. (2009), i.e. an axisymmetric 


Table 1 The list of galaxy component parameters. 


Component 

MIMq] 

a [kpc] 

b [kpc] 

Bulge 

1.4 X 10“ 

0.0 

0.3 

Disk 

9.0 X 10“ 

3.3 

0.3 

Halo 

7.0 X 10“ 

0.0 

25.0 


three-component model, where the bulge, disk, and halo are 
described by Plummer-Kuzmin models (Miyamoto & Nagai 
1975) with the potential 



The parameters a, b, and M of the Milky Way model are 
given in Table 1 for the three components. For details of 
the rotation curve, tidal field and the saddle points of the 
effective potential see Ernst et al. (2010); Just et al. (2009). 

We cover a large range in particle number N and Roche 
volume filling factor A. For the analysis of the mass loss we 
take the average over sets of random realisations in order to 
reduce the impact of random noise. For details see Section 
3. As the main parameter to measure the dissolution time 
we use the half-number-time t^o, where 50% of the initial 
particles are lost. 

For the solution of the A^-body problem the W-body 
programs NBODY6tid and lyS-GRAPEH-GPU were used. 

For this study, we use A', the 99% Roche volume fill¬ 
ing factor, as a measure for the Roche volume filling since 
the 99% Lagrange radius is a statistically more robust mea¬ 
sure than the 100% Lagrange radius. For King (1966) mod¬ 
els, which have a cutoff radius where the density drops to 
zero, the ratios between 99% and 100% Lagrange radius 
are fixed. The conversion factor between rg 9 % and rigo% 
is 1.521 for a ILo = 6 King model. The conversion factor 
between r/j and rgg% is 4.486 for a Wg = 6 King model. 

After the random realization of positions, velocities and 
stellar masses the Jacobi radius rj and the Lagrange radii 
rj% were determined, the cluster size was scaled to realize 
the selected Roche volume filling factors A' in the tidal field 
(Table 2), before the simulation was started. 

2.1 nbody6tid code 

Originally, the program NBODY6tid was written for Galac¬ 
tic centre studies (Ernst, Just & Spurzem 2009; Ernst 
2009) and called NB0DY6GC. Later the three-component 
Plummer-Kuzmin Milky Way model based on Eqn. (2) 
and Table 1 was added to treat the tidal field similar 
to (^-GRAPEH-GPU (see below), and the program was re¬ 
named nbody6tid. In nbody6tid, the galactic centre po¬ 
sition is modelled as a pseudo-particle carrying the Galac¬ 
tic potential in an orbit around the star cluster, which 
must be located close to the origin of coordinates. In ad¬ 
dition to the equations of motion of the star cluster N- 
body problem in the comoving coordinate system, which 
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are solved using a fourth-order Hermite integration scheme 
(Makino & Aarseth 1992) with individual (hierarchical) 
block time steps (Aarseth 2003), NBODY6tid solves the 
equations of motion for the galactic centre orbit with a time- 
transformed eighth-order composition scheme (McLachlan 
1995; Mikkola & Tanikawa 1999; Mikkola & Aarseth 2002; 
Preto & Tremaine 1999; Yoshida 1990). The tidal force of 
the background potential acts on all particles in the A^-body 
system. It is added to the regular force part of NBODY6 
(Aarseth 2003; Ahmad & Cohen 1973). Furthermore, the 
tidal force is added as a perturbation to the KS regulariza¬ 
tion (Kustaanheimo & Stiefel 1965) of NBODY6 (Aarseth 
2003). The time derivative of tidal acceleration (“jerk”) is 
also calculated and added appropriately in the fourth-order 
Hermite integration scheme. Also, the Chandrasekhar dy¬ 
namical friction force (Binney & Tremaine 1987, 2008; 
Chandrasekhar 1943) is implemented using an implicit mid¬ 
point method (Mikkola & Aarseth 2002). Stellar evolution 
is modelled with the fitting formulas of Hurley, Pols & 
Tout (2000). Kicks for compact objects are applied in the 
NBODY6tid runs. For stellar mass black holes and neu¬ 
tron stars the kick velocities are drawn from a Maxwellian 
with a ID dispersion which equals twice the velocity unit 
(Heggie & Mathieu 1986) in km s“^ while for the white 
dwarfs the corresponding ID dispersion is chosen to be 5 
km s“^ (e.g. Fellhauer et al. 2003, for a lower bound). For 
high particle numbers N = 20k, 50k the serial GPU variant 
NBODY6tidgpu is partly used in the present study, based 
on NB0DY6GPU (Nitadori & Aarseth 2012), a version of 
NB0DY6 (Aarseth 2003) which uses NVIDIA type GPU’s 
with CUDA' library support for the force and jerk evalua¬ 
tions. 

2.2 (^-GRAPEH-GPU code 

We also show a few models calculated with the ip- 
GRAPEH-GPU A^-body code which also uses the fourth-order 
Hermite integration scheme (Makino & Aarseth 1992) with 
individual (hierarchical) block time steps and includes the 
three-component Plummer-Kuzmin model based on Eqn. 
(2) and Table 1 . The external gravity part for the tidal field 
is calculated in the galactocentric reference frame. Against 
rounding errors the W-body problem of the star cluster is 
solved in the local cluster frame as in NBODY6tid. The 
first version of the code was written from scratch in C and 
originally designed to use the GRAPE6a clusters for the N- 
body task integration (Harfst et al. 2007). In the present 
version of the (/?-GRAPEH-GPU code we use the NVIDIA 
type GPU’s with CUDA library support, with the exter¬ 
nal SAPPORO library (Gaburov, Harfst & Portegies Zwart 
2009), which emulates for us the standard GRAPE6a library 
calls on the NVIDIA GPU hardware. The (p-GRAPEH-GPU 
code was extensively tested and already long time success¬ 
fully used in our earlier Milky Way star cluster dynamical 
mass loss simulations; Ernst et al. (2010); Just et al. (2009); 

* http://www.nvidia.com 


Table 2 Overview of the parameter space and of the ensem¬ 
bles. Each series comprises a set of ensembles with different parti¬ 
cle numbers N, NBODYbTID ensemble sizes n or yj-GRAPE-l-GPU 
ensemble sizes n' . The quantities Wo and Rg are the 

Roche volume filling factor, the King parameter and the galac¬ 
tocentric radius, respectively. The letter “k” stands for “kilo” 
(=1000). A “-I-” or a means that the corresponding ensemble 
has been calculated with NBODYbxiD or with (p-GRAPE-l-GPU, re¬ 
spectively. A means that the corresponding ensemble has not 
been calculated. The number in round brackets after a denotes 
the number of :p-GRAPE-l-GPU simulations, i.e. the :p-GRAPE-l-GPU 
ensemble size n'. 


Series 

F 

Ui 

U 2 

Oi 

02 

03 

Wo 

6 

6 

6 

6 

6 

6 

A = rh/rj 

0.149 

0.111 

0.074 

0.222 

0.446 

0.669 

A' = rggyjrj 

2/3 

1/2 

1/3 

1 

2 

3 

A = rioo%/rj 

1 

3/4 

1/2 

3/2 

3 

9/2 

Rg [kpc] 

8 

8 

8 

8 

8 

8 

{N,n) 

F 

Ui 

U 2 

Oi 

O 2 

O 3 

(50,512) 

+ 

+ 

+ 

+ 

+ 

+ 

(100,256) 

+ 

+ 

+ 

+ 

+ 

+ 

(200,128) 

+ 

+ 

+ 

+ 

+ 

+ 

(500, 64) 

+ 

+ 

+ 

+ 

+ 

+ 

(Ik, 32) 

+ 

+ 

+ 

+ 

+ 

+ 

(2k, 16) 

+ 

+ 

+ 

+ 

+ 

+ 

(5k, 8) 

+ 

+ 

+ 

+ 

+ 

+ 

(10k, 4) 

+ 

+ 

+ 

+ 

+ 

+ 

(20k, 2) 

+ 

- 

*(3) 

+* (3) 

+* (3) 

+* (6) 

(50k, 1) 

- 

- 

- 

+*(1) 

*(3) 

*(6) 

(100k,-) 

- 

- 

- 

- 

*(6) 

*(6) 

(1000k, -) 

- 

- 

- 

*(1) 

- 

- 


Kharchenko et al. (2009).^ The stellar evolution treatment 
of (/5-GRAPEH-GPU is in detail described in Kharchenko et 
al. (2009). Eor the (/)-GRAPEH-GPU runs in the present study 
kicks for compact objects have not been applied for techni¬ 
cal reasons. 

3 Parameter space 

Table 2 shows an overview over the parameter space. We 
employ ensemble averaging over n NBODY6tid simula¬ 
tions or n' (/3-GRAPEH-GPU simulations. 

We have chosen the minimum particle number to be 
N = 50. We remark that the definition of the two-body 
relaxation time in Eqn. (3) below breaks down for small 
N due to the minimum in N/ ln( 7 A^) at W = 136 for 
7 = 0.02. The evolution of the system is no longer gov¬ 
erned by small-angle scatterings in such a small-W regime. 

All models were King (1966) models with the King pa¬ 
rameter Wq = 6 placed on a circular orbit at the Galacto¬ 
centric radius Rg = 8 kpc in 2 ; = 0 plane of the Galac¬ 
tic tidal field based on Eqn. (2) and Table 1. In all mod- 

^ The first original public version of the ip-GRAPE+GPU code can found 
here: 

ftp://ftp.mao.kiev.ua/pub/users/berczik/phi-GRAPE+GPU/ 
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els, we applied a Kroupa (2001) initial mass function with 
0.08 < tti/Mq < 100.00, where m is the stellar mass. 
We remark that the IMF can have a drastic effect on the 
lifetimes of star clusters as Engle (1999) demonstrated nu¬ 
merically. For the metallicity we chose the solar metallicity 
Z = 0.02 in all models. We remark that this value may not 
be up-to-date anymore. However, it allows for comparisons 
with older simulations. 

The lower part of Table 2 shows an overview over the 
simulations which have been carried out. A “H-” or a 
means that the corresponding ensemble has been calculated 
with nbody6tid or with (^-grapeh-gpu, respectively. The 
ensembles marked with a sign have not been simulated. 
The number in round brackets after a denotes the num¬ 
ber of additional (^-GRAPEH-GPU simulations, i.e. the (p- 
GRAPEH-GPU ensemble size n'. 

We remark that Baumgardt & Makino (2003) investi¬ 
gated the Roche volume hlling case with A = 1 correspond¬ 
ing to our series F. 

3.1 nbody6tid parameter space 

The following parameter space was adopted for the 
NBODY6tid simulations: For the ordered pairs {N,n) 
we choose (50,512), (100,256), (200,128), (500,64), 
(1000,32), (2000,16), (5000,8), (10000,4), (20000,2), 
(50000,1), where N is the particle number and n the 
number of NBODY6tid runs per ensemble. For each or¬ 
dered pair {N, n) we computed 6 ensembles correspond¬ 
ing to the Roche volume hlling factors A' = rgg%lrj = 
1/3 (172), 1/2 (C/i),2/3 (F),l (Oi),2 (02),3 (Og)-The 
full NBODY6tid parameter space comprised approx. 50 en¬ 
sembles or approx. 6000 runs, respectively. 

3.2 (/5-grapeh-gpu parameter space 

The (p-GRAPEH-GPU simulations were computed addition¬ 
ally. The main reason to use both (similar) W-body codes 
was, that we wanted to compare the results of both programs 
and to apply the advantages of both programs to the same 
star cluster dynamical evolution scientihc application. The 
(/5-GRAPEH-GPU simulations include one model with = 1 
million for series Oi. 

4 Results 

4.1 Random scatter 

For a cluster with N = 1000 or smaller there is a consid¬ 
erable scatter in the realized cluster mass and size due to 
the random population of the high mass end of the IMF and 
of the outer shells of the cluster. Figure 1 shows the time 
evolution of all 128 runs of the ensemble with N = 200 of 
Series Oi. The ensemble mean is marked by the thick solid 
black line. There is a strong scatter in mass loss times. Fig¬ 
ure 2 shows the half-number times fso as function of initial 



0 500 1000 1500 2000 

t [Myr] 


Fig. 1 Example for the evolution of the particle number, ensem¬ 
ble with (A, n) = (200,128). The ensemble mean is marked by 
the thick solid black line. 



0 12 3 4 5 6 

log(Mo [Ms„|]) 


Fig. 2 Half-number times versus initial mass for all runs of Se¬ 
ries Oi. The particle numbers are marked by labels. The mean 
values are marked by diamonds. 

mass for all (A, n) ensembles of series 0\. For each hxed 
N there is a clear anticorrelation with initial mass showing 
that a few high-mass stars accelerate the dissolution signih- 
cantly. Figure 3 shows the corresponding distributions of 
for the low-A ensembles of series Oi. The binsize is 0.02 
dex in log]^Q(f 5 o). The ensemble median Q^q is marked by 
the solid black line. The dotted and dashed lines mark the 
30% and 70% quantiles Q 30 and Q^o of the corresponding 
distribution. The COYOTE library in IDL (Fanning 2011) has 
been used. 

In Figure 4 we compare the evolution of the half-mass 
radius as function of relative mass loss M(t) /Mq in terms 
of A for different series. This shape parameter A increases 
during the evolution similar to the result of Fukushige & 
Heggie (1995) and it is independent of N. There is a clear 
separation of the different series showing that there is a 
memory of the initial relative size of the half-mass radius. 

4.2 Mass loss 

The left-hand side of Figure 5 shows the time evolution 
of the arithmetic ensemble mean of the normalized parti- 


www.an-journal.org 


© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 





794 


A. Ernst, P. Berczik, A. Just, T. Noel: Dissolution of open star clusters 


30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 
30 

z 20 
10 
0 



Series 0^, N = 1 00 



Series 0^, N = 200 



Series 0-,, N = 500 



Series 0^, N = 1 000, 



Series 0 -,, N = 2000 



Series 0^, N = 5000 


Series 0^, N = 1 0000 



2.0 2.5 3.0 

°9lo(t5o) 


Fig. 3 Distributions of half-number times tso for the low-A^ en¬ 
sembles of series Oi. The binsize is 0.02 dex in log]^g(f 5 o). The 
ensemble median Qso is marked by the solid black line. The dot¬ 
ted and dashed lines mark the 30% and 70% quantiles Q 30 and 
Q 70 of the corresponding distribution, respectively. 



M/Mo 


Fig. 4 Evolution of the half-mass Roche volume filling factor A 
as a function of the bound mass M/Mq. 


cle number N{t)/NQ within three times the Jacobi radius, 
which we defined to be bound. 

All curves decline monotonically as the simulated star 
clusters lose mass. The moderately Roche volume underfill¬ 
ing series Ui and the strongly Roche volume underfilling se¬ 
ries U 2 dissolve slower than the Roche volume filling series 
F. The reason is that they first need a phase of expansion 
to fill the Roche volume (e.g. Engle 1999). It can be seen 
that in the moderately Roche volume overfilling series O 2 
the iV-dependence is much weaker than in the Roche vol¬ 
ume overfilling series Oi. Moreover, in the strongly Roche 
volume overfilling series O 3 the fV-dependence has almost 
completely vanished. This indicates that two-body relax¬ 
ation is not responsible for the dissolution. 

The right-hand side of Figure 5 shows the time evolution 
of the arithmetic ensemble mean of the normalized mass 
M(f)/Mo within three times the Jacobi radius. All curves 
show a strong initial decrease due to the stellar evolution 
mass loss. The half-mass times are significantly shorter than 
the half-number times due to the stellar evolution mass loss. 
As a consequence the binding energy and Roche volume 
decrease, which depends on the cluster mass, faster than 
2 -body relaxation, which depends on the number of stars. 
These differences decrease with increasing Roche volume 
filling factor. 

4.3 Dissolution times 

There are four relevant time scales involved in the dissolu¬ 
tion of star clusters in the Galactic tidal field: (i) The two- 
body relaxation time fj-x, (ii) crossing time ter, (iii) orbital 
time forb and (iv) stellar evolution time fgtev The first three 
of them scale as 

N 


icrOC 


(4) 
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Fig. 5 Evolution of the particle number within three times the Jacobi radius (left panels) and cluster mass within three times the Jacobi 
radius (right panels). From left to right in each panel: N = 200, 500, Ik, 2k, 5k, lOfc, 20k, 50k. The sharp initial decrease in the total 
mass is due to the stellar evolution mass loss. 
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Tables The slopes 0 : 50 . 


Ser. 

7 

F 

C/i 

C/2 

Oi 

O 2 

O 3 

*50 

0.02 

0.708 

0.843 

0.872 

0.573 

0.238 

0.100 

*50 

0.11 

0.630 

0.752 

0.764 

0.494 

0.212 

0.100 


^orb 


GMg 


-1/2 


(5) 


where Mg , Rg , and 7, are the the enclosed Galaxy mass, the 
galactocentric radius, and the factor in the Coulomb log¬ 
arithm (Giersz & Heggie 1994, 1996), respectively, rj is 
given by Eqn. (1). We note that we use the initial two-body 
relaxation time Rx throughout this study and not the current 
one. The stellar evolution time Gtev depends only on the 
IMF and the metalllicity, which are fixed in our study. For 
ensembles with large particle numbers N and correspond¬ 
ingly large two-body relaxation times Rx the mass loss due 
to stellar evolution is clearly distinguishable from the two- 
body relaxation driven evolution and becomes important 
with respect to it as can be seen in right-hand side panels 
of Figure 5. 


Figure 6 shows that, in the low-JV regime of OCs, the 
half-number time scales directly with a power (7, A') 
of the two-body relaxation time Gx- Figure 6 also shows 
the corresponding straight-line fits for the determination of 
2^50(7) A'). The 6 lowest-data points of the NBODY6tid 
ensembles have been used for the fitting of power laws. 
For the least-squares-fitting, we used the MPFIT package 
in IDL (Markwardt 2009; More 1978, for the Fevenberg- 
Marquardt algorithm). 

We find in this study for the half-number time t^o the 
approximate expression 

^50 

T 

with T « 125 Myr cx; forb and G = 80 — 100. The expo¬ 
nents 


1 N 
C ln(77V) 


(6) 


d log fso 
2-50 ,, , 

d log trx 


(7) 


for the half-number times are given in the legend of Figure 
6 and Table 3, and they are plotted in Figure 7 against the 
Roche volume filling factor. 

The upper and lower errors in Figures 6 , 8 and 9 are 
given by 


A+ logio ti,N = logio(Q70,i,Af) — logio(Q50,i,Af)5(8) 
^ loglO^i.Af = logio(Q50,i,Af) ~ loglo(Q30,i.Af)5(9) 

where Q 3 o,n, Qso.Af and Qyq^n are 30%, 50% and 70% 
quantiles of the corresponding distribution of dissolution 
times ti. The quantiles have been calculated with an IDL 


1,0 
0.8 

o 

0.4 
0.8 
0.0 

0 12 3 4 

^ “^99%/^’j 

Fig. 7 Exponent x^o from Eqn. (6) as a function of A' = 
rgg%/rj. Thick dots: 7 = 0.02; thin dots: 7 = 0.11. 


7 = 0.02 
7 = 0.11 


routine by Hong, Schlegel & Grindlay (2004). The weights 
used in the fitting procedure are given by 


= 1 /[|A+logiofi,Ar| -f |A logigfi^jvl]^ • (10) 

The dependence of 0 : 50 ( 7 , A') on the 99% Roche vol¬ 
ume filling factor A' and the 7 factor in the Coulomb log¬ 
arithm can be seen in Figure 7. We calculated the scaling 
exponents for two different values of the 7 parameter in 
the Coulomb logarithm: 7 = 0.02 (Giersz & Heggie 1996, 
multi-mass case) and 7 = 0.11 (Giersz & Heggie 1994, 
equal-mass case only for comparison) to show the differ¬ 
ence between these two cases. 

Figure 8 shows example fits for the times ^ 20 , ^so (half¬ 
number time) and tso for Series F", C/ 2 , Oi, O 2 . Here we 
used only 7 = 0.02 in the Coulomb logarithm. For the least- 
squares-fitting, we used the MPFIT package in IDL (Mark¬ 
wardt 2009; More 1978, for the Fevenberg-Marquardt algo¬ 
rithm). 

Figure 9 shows all scaling exponents 


d\og{ti) 

filog(Gx) 


( 11 ) 


for i = 1 — 90 (in percent), together with the correspond¬ 
ing values. The time ti is defined as the time when the 
cluster has lost i% of its initial particle number. We calcu¬ 
lated the scaling exponents for two different values of the 7 
parameter in the Coulomb logarithm. We used 


/2 

\red,i 


1 

Ndol 


Ndoi 






( 12 ) 


with the weights of Fqn. (10), where is the value of 
the fitted power law function Eqn. ( 6 ), the yi^ are the data 
and Adot is the number of degrees of freedom (Markwardt 
2009). From Figure 9 we can see that the power law index 
Xi is a weak function of the mass loss fraction. The x'red i 1® 
not easy to interpret since the quantiles of the correspond¬ 
ing statistic with the weights of Eqn. (10) are not known. 
However, one can gain an insight regarding the relative be¬ 
haviour of x'red i 1°'' different fi’s. It can be seen that is 
a robust measure for the dissolution time. 
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Fig. 6 Half-number times in Myr as functions of N/ ln(0.02Af) (dots) and straight-line fits, each using the six lowest-A^ data points 
(for N = 200, 500, Ik, 2k, 5k, 10k) to obtain the slopes in the low-A" OC limit. The results of a few g>-GRAPE-l-GPU simulations are 
also shown. The gi-GRAPE-l-GPU simulation with 1 million particles belongs to series 0\. The weights and error bars are calculated from 
the quantiles Qsq^n, Qso.iv and Q 7 o,n (see text and Eqns. (8) - (10)). 


Table 4 Averaged parameters (logjg(f 5 o [Myr]),«:) from fits 
with Equation (15). 


Series 

logio(f50 [Myr]) 

K 

O 2 

2.30 ±0.09 

2.25 ±0.28 

O 3 

2.16 ±0.04 

2.18 ±0.20 


4.4 Log-logistic growth 


d\nNe / NA N 

d\nt No) '^No' 

The solution is given by 


(13) 


jVe(t) ^ 

No tlo + 1- 


(14) 


The evolution of N {t)/No can then be described by the law 


Figure 10 shows the time evolution of the escaper fraction 
Ne(t)/No = 1 — N(t)/No with a logarithmic time axis, 
where N{t) (and Nq) are taken to be the current (initial) 
particle number within three times the Jacobi radius. 

In the Roche volume overfilling limit the evolution of 
the escaper fraction can, at least empirically, be approxi¬ 
mately described by a log-logistic differential equation in 
logarithmic time, 


‘'50 

t^o + t- 


(15) 


Njt) f? 

No 1-50 

The best-fit exponent k and the best-fit half-number time 
(50 can be determined for the parameter space covered in 
this study. Table 4 shows the parameters of least-quares fits 
with the Equation 


logf = logfso -b - log 

K 


Ne/Np 
1 - Ne/No 


(16) 
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Fig.8 Example fits. From left to right: For f 2 o, tso, fso- From top to bottom: For Series F, U 2 , Oi, O 2 ■ The fits of t 2 o can be biased 
due to stellar evolution mass loss. We used 7 = 0.02 in the Coulomb logarithm. The weights and error bars are calculated from the 
quantiles Qao.iv, Qso.iv and Q 7 o,n (see text and Eqns. ( 8 ) - (10)). 


For the least-squares-fitting, we used the MPFIT package 
in IDL (Markwardt 2009; More 1978, for the Levenberg- 
Marquardt algorithm). The fits are shown in the two lowest 
panels on the right-hand side of Figure 10 (for series O 2 and 
O 3 ). For the series shown in the upper panels we suspect a 
transition from log-logistic to logistic with some other con¬ 
tribution from left to right, where the factor 1 Jt is due to the 
stellar evolution (Earners, Baumgardt & Gieles 2010) and 
the other contribution is due to the relaxation-driven evolu¬ 
tion. 

There are two competing processes for populating the 
potential escaper reservoir above the critical Jacobi energy. 
Firstly, scattering by stellar encounters scales with the re¬ 
laxation time and depends on the particle number N. Sec¬ 
ondly, cluster mass loss by stellar evolution and by escaping 
stars lowers the cluster potential well and lifts the criticalJa- 
cobi energy, which shifts new stars above the critical value. 
The fractional mass loss rate is independent of N and re¬ 
sults in a particle loss rate proportional to N, to the num¬ 
ber of escapers N,. and via stellar evolution to 1 /t (Tamers, 
Baumgardt & Gieles 2010). For the overfilling clusters (se¬ 


ries O 2 , O 3 ) two-body encounters are negligible leading to 
the log-logistic behaviour. For the more concentrated clus¬ 
ters, where two-body encounters play an important role, the 
interplay of the different timescales is discussed in detail in 
Tamers, Baumgardt & Gieles (2010). Only for very large 
N, where the relaxation time is very long, the factor 1 /t by 
stellar evolution in the dissolution timescale becomes rel¬ 
evant again. This leads to a reduction of t^o with respect 
to the power law dependence derived in Eqn. ( 6 ) for small 
filling factors. 

5 Conclusions 

We have carried out a parameter study of open star clusters 
with the parameters {N,X = rtjrj). We have found the 
following results: 

1. The A^-dependence of the dissolution time in units of the 
two-body relaxation time is well fitted by a power law. 
The power law index is a function of the Roche volume 
filling factor A and the 7 factor in the Coulomb loga¬ 
rithm. It decreases with increasing A with the limiting 
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Fig. 9 Top panel: Scaling exponents Xi = dlog(ti)/dlog(trh)- The time ti is defined as the time when the cluster has lost i% of its 
initial particle number. We used a median smoothing with a smoothing width of 11 in i. Bottom panel: The corresponding x'red,i values. 
For low i we find large x'red i due to the stellar evolution mass loss. The half-number time seems to be a robust measure with respect to 
the x?ed i values. Thick lines: 7 = 0.02; thin lines: 7 = 0.11. The weights are calculated from the quantiles Qao.jv, Qso.iv and Qro.iv 
(see text and Eqn. (10)). 


value of zero in the overfilling limit. Particularly in the 
underfilling limit, the power law index has been found 
to depend on the value of the 7 parameter adopted in the 
Coulomb logarithm. 

2. Our study suggests that open star clusters in the Roche 
volume overfilling regime dissolve mainly due to the 
changing cluster potential by lifting stars above the de¬ 
creasing critical Jacobi energy. We call this mechanism 
“mass-loss driven dissolution” in contrast to the “two- 
body relaxation driven dissolution” which occurs from 
the Roche volume underfilling regime up to the Roche 
volume filling case (see also Whitehead et al. 2013, 
based on simpler models). 

3. In the Roche volume overfilling limit the escaper frac¬ 
tion Ne{t)/No obeys, at least empirically, approxi¬ 
mately a log-logistic differential equation in logarithmic 
time. 

We make the following remarks: 

1. The mass-loss driven dissolution provides a mechanism, 
which is responsible for the dissolution of OCs which 
have survived the gas expulsion phase with a relatively 
large initial half-mass radius as observed for example in 


the Pleiades cluster (Converse & Stabler 2010). It is a vi¬ 
able mechanism besides the dissolution due to encoun¬ 
ters with giant molecular clouds (Wielen 1971, 1985). 

2. The mass-loss driven dissolution in combination with a 
large scatter in Roche volume filling factors (Ernst & 
Just 2013) also naturally explains the large scatter in 
the lifetimes of open clusters (Wielen 1971). In a fu¬ 
ture paper we plan to investigate the impact of the newly 
found A^-independence of the dissolution timescale for 
the overfilling clusters on the CFR using the observed 
mass-age distribution of OCs in the solar neighbour¬ 
hood. 

3. Due to the intricacy of the problem a detailed theory of 
mass-loss driven dissolution, which explains the transi¬ 
tion from the over- to underfilling limits in terms of the 
dependence of 0:50 on A as shown in Figure 7 has not yet 
been developed and is beyond the scope of the present 
experimental study. 
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A Resonance condition 


It is possible to write down a resonance condition (Ernst & Just 
2013), 


m _ _ / A- f Ges\ 

n Ores V 2 ) \ rj J 


(Al) 



Fig. B1 Stability curve for a Plummer model and King models. 
Value of rh /ry for Wo = 11 by Giirkan and Rreitag, priv. comm. 


where m and n are natural numbers and the orbital frequencies 
^res ~ ^('^res ) and G = 2 'k are related to the orbital time 
of a star in the star cluster and the orbital time of the star cluster 
on a circular orbit around the galaxy, respectively. Pc is the ratio 
between epicyclic and circular frequency and r^es is a resonance 
radius. For a flat rotation curve we have Pc = \/2- For the Milky 
Way model with the parameters given in Table 1 we have Pc ~ 
1.37 at i?g = 8 kpc. Therefore at certain values of A = ru/rj 
resonance effects may occur and play a role in the evolution. 

The central periodic orbit in the largest regular (i.e. non- 
chaotic) island in the Poincare surfaces of section at the critical 
Jacobi energy has (rres/fj) = 0.345. This island corresponds 
to quasiperiodic retrograde orbits (see also Ernst et al. 2008; 
Fukushige & Heggie 2000). It may be that the corresponding res¬ 
onance is linked with the properties of star clusters and connected 
with the occurence of two discrete types of star clusters, open and 
globular clusters. Moreover, the location of the resonance may be 
used to calibrate the relation (Al) more precisely for star clusters, 
i.e. to determine rres as a function of rn- 


B Stability curve 


In connection with the Roche volume filling factor a stability curve 
can be derived for star clusters. Such a curve has been derived 
for the first time in Fukushige & Heggie (1995, their Eqn. 28; see 
also their Figure 13) and may therefore be called the “Fukushige- 
Heggie stability curve” for star clusters. The Jacobi energy of a 
typical star at the half-mass radius is given by 


ej,h 


K W 

Mel Mel 


1GM^3 

2 ^^ 


(Bl) 


with A = rhjrj, where K and W are the kinetic and potential 
energies, respectively. The critical Jacobi energy is given by 


3 GMel 

^J.crit — ^ 

2 rj 

We obtain for a Plummer model with 

IR _ 37r GMl 

2 “ 64 rpi ’ 

and rh/rp\ « 1.305 the stability curve 

t^J,h t^J^erit ^ 1.305 * TT A 
ej.erit 32A 3 

For a King model with Wq = 6 we have 


(B2) 


(B3) 


(B4) 
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K = = 


GMl 

Arv 


(B5) 


with rh/rv ~ 0.804 (Half-mass radius in Ai-body units G = 
Mci = —AE = 1, Giirkan, Freitag & Rasio 2004, Table 1), where 
E is the total energy. We obtain the stability curve 


ej,crit 0.804 A /D/T'. 

--f ^— 4 (Bbj 

6j,crit OA O 

The stability curves in Eqns. (B4) and (B6) are shown in Figure 
B1 for a Plummer model and 3 King models and may be also con¬ 
nected with the occurence of two discrete types of star clusters, 
open and globular clusters. The value of this function measures 
the relative Jacobi energy difference between a typical star at the 
half-mass radius and the critical Jacobi energy. Note that Eqn. (B2) 
and the virial relations in Eqns. (B3) and (B5) break down for large 
A, i.e. only the first zero at small A may be of physical significance. 
We expect relaxation driven dissolution in the regime to the left of 
the first zero and mass-loss driven dissolution in the regime to the 
right of the first zero. 
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